home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlatrd.f < prev    next >
Text File  |  1997-06-25  |  10KB  |  260 lines
  1.       SUBROUTINE DLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     October 31, 1992
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          UPLO
  10.       INTEGER            LDA, LDW, N, NB
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), E( * ), TAU( * ), W( LDW, * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DLATRD reduces NB rows and columns of a real symmetric matrix A to
  20. *  symmetric tridiagonal form by an orthogonal similarity
  21. *  transformation Q' * A * Q, and returns the matrices V and W which are
  22. *  needed to apply the transformation to the unreduced part of A.
  23. *
  24. *  If UPLO = 'U', DLATRD reduces the last NB rows and columns of a
  25. *  matrix, of which the upper triangle is supplied;
  26. *  if UPLO = 'L', DLATRD reduces the first NB rows and columns of a
  27. *  matrix, of which the lower triangle is supplied.
  28. *
  29. *  This is an auxiliary routine called by DSYTRD.
  30. *
  31. *  Arguments
  32. *  =========
  33. *
  34. *  UPLO    (input) CHARACTER
  35. *          Specifies whether the upper or lower triangular part of the
  36. *          symmetric matrix A is stored:
  37. *          = 'U': Upper triangular
  38. *          = 'L': Lower triangular
  39. *
  40. *  N       (input) INTEGER
  41. *          The order of the matrix A.
  42. *
  43. *  NB      (input) INTEGER
  44. *          The number of rows and columns to be reduced.
  45. *
  46. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  47. *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  48. *          n-by-n upper triangular part of A contains the upper
  49. *          triangular part of the matrix A, and the strictly lower
  50. *          triangular part of A is not referenced.  If UPLO = 'L', the
  51. *          leading n-by-n lower triangular part of A contains the lower
  52. *          triangular part of the matrix A, and the strictly upper
  53. *          triangular part of A is not referenced.
  54. *          On exit:
  55. *          if UPLO = 'U', the last NB columns have been reduced to
  56. *            tridiagonal form, with the diagonal elements overwriting
  57. *            the diagonal elements of A; the elements above the diagonal
  58. *            with the array TAU, represent the orthogonal matrix Q as a
  59. *            product of elementary reflectors;
  60. *          if UPLO = 'L', the first NB columns have been reduced to
  61. *            tridiagonal form, with the diagonal elements overwriting
  62. *            the diagonal elements of A; the elements below the diagonal
  63. *            with the array TAU, represent the  orthogonal matrix Q as a
  64. *            product of elementary reflectors.
  65. *          See Further Details.
  66. *
  67. *  LDA     (input) INTEGER
  68. *          The leading dimension of the array A.  LDA >= (1,N).
  69. *
  70. *  E       (output) DOUBLE PRECISION array, dimension (N-1)
  71. *          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
  72. *          elements of the last NB columns of the reduced matrix;
  73. *          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
  74. *          the first NB columns of the reduced matrix.
  75. *
  76. *  TAU     (output) DOUBLE PRECISION array, dimension (N-1)
  77. *          The scalar factors of the elementary reflectors, stored in
  78. *          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
  79. *          See Further Details.
  80. *
  81. *  W       (output) DOUBLE PRECISION array, dimension (LDW,NB)
  82. *          The n-by-nb matrix W required to update the unreduced part
  83. *          of A.
  84. *
  85. *  LDW     (input) INTEGER
  86. *          The leading dimension of the array W. LDW >= max(1,N).
  87. *
  88. *  Further Details
  89. *  ===============
  90. *
  91. *  If UPLO = 'U', the matrix Q is represented as a product of elementary
  92. *  reflectors
  93. *
  94. *     Q = H(n) H(n-1) . . . H(n-nb+1).
  95. *
  96. *  Each H(i) has the form
  97. *
  98. *     H(i) = I - tau * v * v'
  99. *
  100. *  where tau is a real scalar, and v is a real vector with
  101. *  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
  102. *  and tau in TAU(i-1).
  103. *
  104. *  If UPLO = 'L', the matrix Q is represented as a product of elementary
  105. *  reflectors
  106. *
  107. *     Q = H(1) H(2) . . . H(nb).
  108. *
  109. *  Each H(i) has the form
  110. *
  111. *     H(i) = I - tau * v * v'
  112. *
  113. *  where tau is a real scalar, and v is a real vector with
  114. *  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
  115. *  and tau in TAU(i).
  116. *
  117. *  The elements of the vectors v together form the n-by-nb matrix V
  118. *  which is needed, with W, to apply the transformation to the unreduced
  119. *  part of the matrix, using a symmetric rank-2k update of the form:
  120. *  A := A - V*W' - W*V'.
  121. *
  122. *  The contents of A on exit are illustrated by the following examples
  123. *  with n = 5 and nb = 2:
  124. *
  125. *  if UPLO = 'U':                       if UPLO = 'L':
  126. *
  127. *    (  a   a   a   v4  v5 )              (  d                  )
  128. *    (      a   a   v4  v5 )              (  1   d              )
  129. *    (          a   1   v5 )              (  v1  1   a          )
  130. *    (              d   1  )              (  v1  v2  a   a      )
  131. *    (                  d  )              (  v1  v2  a   a   a  )
  132. *
  133. *  where d denotes a diagonal element of the reduced matrix, a denotes
  134. *  an element of the original matrix that is unchanged, and vi denotes
  135. *  an element of the vector defining H(i).
  136. *
  137. *  =====================================================================
  138. *
  139. *     .. Parameters ..
  140.       DOUBLE PRECISION   ZERO, ONE, HALF
  141.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 )
  142. *     ..
  143. *     .. Local Scalars ..
  144.       INTEGER            I, IW
  145.       DOUBLE PRECISION   ALPHA
  146. *     ..
  147. *     .. External Subroutines ..
  148.       EXTERNAL           DAXPY, DGEMV, DLARFG, DSCAL, DSYMV
  149. *     ..
  150. *     .. External Functions ..
  151.       LOGICAL            LSAME
  152.       DOUBLE PRECISION   DDOT
  153.       EXTERNAL           LSAME, DDOT
  154. *     ..
  155. *     .. Intrinsic Functions ..
  156.       INTRINSIC          MIN
  157. *     ..
  158. *     .. Executable Statements ..
  159. *
  160. *     Quick return if possible
  161. *
  162.       IF( N.LE.0 )
  163.      $   RETURN
  164. *
  165.       IF( LSAME( UPLO, 'U' ) ) THEN
  166. *
  167. *        Reduce last NB columns of upper triangle
  168. *
  169.          DO 10 I = N, N - NB + 1, -1
  170.             IW = I - N + NB
  171.             IF( I.LT.N ) THEN
  172. *
  173. *              Update A(1:i,i)
  174. *
  175.                CALL DGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
  176.      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
  177.                CALL DGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
  178.      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
  179.             END IF
  180.             IF( I.GT.1 ) THEN
  181. *
  182. *              Generate elementary reflector H(i) to annihilate
  183. *              A(1:i-2,i)
  184. *
  185.                CALL DLARFG( I-1, A( I-1, I ), A( 1, I ), 1, TAU( I-1 ) )
  186.                E( I-1 ) = A( I-1, I )
  187.                A( I-1, I ) = ONE
  188. *
  189. *              Compute W(1:i-1,i)
  190. *
  191.                CALL DSYMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
  192.      $                     ZERO, W( 1, IW ), 1 )
  193.                IF( I.LT.N ) THEN
  194.                   CALL DGEMV( 'Transpose', I-1, N-I, ONE, W( 1, IW+1 ),
  195.      $                        LDW, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
  196.                   CALL DGEMV( 'No transpose', I-1, N-I, -ONE,
  197.      $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
  198.      $                        W( 1, IW ), 1 )
  199.                   CALL DGEMV( 'Transpose', I-1, N-I, ONE, A( 1, I+1 ),
  200.      $                        LDA, A( 1, I ), 1, ZERO, W( I+1, IW ), 1 )
  201.                   CALL DGEMV( 'No transpose', I-1, N-I, -ONE,
  202.      $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
  203.      $                        W( 1, IW ), 1 )
  204.                END IF
  205.                CALL DSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
  206.                ALPHA = -HALF*TAU( I-1 )*DDOT( I-1, W( 1, IW ), 1,
  207.      $                 A( 1, I ), 1 )
  208.                CALL DAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
  209.             END IF
  210. *
  211.    10    CONTINUE
  212.       ELSE
  213. *
  214. *        Reduce first NB columns of lower triangle
  215. *
  216.          DO 20 I = 1, NB
  217. *
  218. *           Update A(i:n,i)
  219. *
  220.             CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
  221.      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
  222.             CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
  223.      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
  224.             IF( I.LT.N ) THEN
  225. *
  226. *              Generate elementary reflector H(i) to annihilate
  227. *              A(i+2:n,i)
  228. *
  229.                CALL DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1,
  230.      $                      TAU( I ) )
  231.                E( I ) = A( I+1, I )
  232.                A( I+1, I ) = ONE
  233. *
  234. *              Compute W(i+1:n,i)
  235. *
  236.                CALL DSYMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
  237.      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
  238.                CALL DGEMV( 'Transpose', N-I, I-1, ONE, W( I+1, 1 ), LDW,
  239.      $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
  240.                CALL DGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
  241.      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
  242.                CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA,
  243.      $                     A( I+1, I ), 1, ZERO, W( 1, I ), 1 )
  244.                CALL DGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
  245.      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
  246.                CALL DSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
  247.                ALPHA = -HALF*TAU( I )*DDOT( N-I, W( I+1, I ), 1,
  248.      $                 A( I+1, I ), 1 )
  249.                CALL DAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
  250.             END IF
  251. *
  252.    20    CONTINUE
  253.       END IF
  254. *
  255.       RETURN
  256. *
  257. *     End of DLATRD
  258. *
  259.       END
  260.